Ultracold bosons in disordered superlattices: Mott- insulators induced by tunneling 
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We analyse the phase diagram of ultra-cold bosons in a one-dimensional superlattice potential with disorder 
using the time evolving block decimation algorithm for infinite sized systems (iTEBD). For degenerate potential 
energies within the unit cell of the superlattice loophole-shaped insulating phases with non-integer filling emerge 
with a particle-hole gap proportional to the boson hopping. Adding a small amount of disorder destroys this 
gap. For not too large disorder the loophole Mott regions detach from the axis of vanishing hopping giving rise 
to insulating islands. Thus the system shows a transition from a compressible Bose-glass to a Mott-insulating 
phase with increasing hopping amplitude. We present a straight forward effective model for the dynamics within 
a unit cell which provides a simple explanation for the emergence of Mott-insulating islands. In particular it 
gives rather accurate predictions for the inner critical point of the Bose-glass to Mott-insulator transition. 

PACS numbers: 



I. INTRODUCTION 

Ultra-cold atomic gases in light induced periodic poten- 
tials have become an important experimental testing ground 
for concepts of solid-state and many-body physics since they 
allow the realization of precisely controllable model Hamil- 
tonians with widely tunable parameters. This development 
was triggered by the theoretical proposal of Jaksch et al. JH] 
that ultra-cold bosonic atoms in an optical lattice are well de- 
scribed by the Bose-Hubbard model and the subsequent obser- 
vation of the superfluid-Mott insulator transition in that sys- 
tem by Greiner et al. J2]. A characteristic feature of light 
induced periodic potentials is the possibility to modify their 
properties in a simple way. E.g. when phase locked lasers 
with different but commensurable frequencies are superim- 
posed to create a periodic dipole potential, different types of 
superlattices with more complex unit cells can be constructed 

OIHIIIHIHtIHI!]. Su P erim P° sed °P tical lattices with non- 
commensurate frequencies furthermore mimic a disordered 
potential ifioll . 

In the present paper we study the phase diagram of ultra- 
cold bosons in a one-dimensional superlattice potential with 
degenerate potential energies and/or degenerate tunneling 
rates within the unit cell. For such a system loophole shaped 
Mott-insulator domains with fractional filling have been pre- 
dicted by Buonsante, Penna and Vezzani f 70 within a multiple- 
site mean-field approach, as well as with a cell strong- 
coupling perturbation approach J^]. In contrast to the Mott 
lobes at integer filling known from the simple Bose-Hubbard 
model, which exist also in the superlattice, the characteristic 
feature of the loophole insulators is a particle-hole gap that 
vanishes at zero boson hopping J. So in the \i — J phase dia- 
gram, where fi is the chemical potential, these domains touch 
the J = line only in a single point. We here perform nu- 
merical simulations using the time evolving block decimation 
algorithm (TEBD) introduced by Vidal et al. ifTTll in the infi- 
nite system variant Il2tl as well as density matrix renormaliza- 
tion group (DMRG) calculations lfl3ll to determine the bound- 
aries of the different Mott-insulating regions in the phase dia- 
gram. We also present a simple effective model that provides a 
straight forward explanation for the emergence of the loophole 



insulators by taking into account hopping processes between 
the sites of degenerate potential energy within a unit cell but 
neglecting tunneling between different unit cells. 

We then study the influence of some additional disorder po- 
tential with continuous, bounded distribution. If the maxi- 
mum amplitude of the disorder is not too large, the Mott lobes 
shrink in a similar way as for the simple Bose-Hubbard model 
ITill . Also the loophole Mott domains shrink. As a con- 
sequence near the critical (fractional) filling Mott-insulating 
islands emerge surrounded by a Bose-glass phase. A rather 
peculiar property of this system is the phase transition from 
a compressible (Bose-glass) phase for small values of the 
bosonic hopping J to an incompressible Mott phase for larger 
tunneling rates, i.e. we have a tunneling induced Mott insula- 
tor. The effective model describing the full dynamics within 
a unit cell provides a simple explanation for this and gives 
good quantitative predictions for the critical value J c for the 
compressible-to-Mott transition. The analytical predictions 
are compared to numerical simulations again using the iTEBD 
algorithm for a superlattice with disorder. 



H. THE MODEL 

We consider ultracold bosonic atoms in an optical superlat- 
tice, having a periodic structure with a period I of some lattice 
sites. As shown in |Q]], the physics of these atoms can be de- 
scribed by the so called superlattice Bose-Hubbard model, ex- 
tensively studied in JH Hfl . We here work mostly in the grand 
canonical BHM, only the calculations of the shape of the loop- 
hole insulators in sections IIIIBI and |IV A| will be performed 
for fixed particle numbers. In second quantisation, the BHM 
reads as 



j 

(1) 



where &j and at are the annihilation and creation operators of 
the bosons at lattice site j, and hj = ajfij is the corresponding 
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number operator. The particles can tunnel from one lattice site 
to a neighbouring one with hopping rate J, tj accounts for the 
variation of the hopping due to the superlattice potential. Vj 
accounts for local variations of the potential energy within a 
unit cell, and fi is the (global) chemical potential. 

A particularly interesting situation arises if there is a degen- 
eracy in the tunneling amplitudes and local potentials within 
the unit cell. This will be studied in detail in the follow- 
ing. For simplicity we focus on a special superlattice structure 
in which only the local potential is varied with period three, 
namely 

v = {v 1 ,v 1 ,v 2 } (2) 

with vi — U< v 2 < Vi < U and t = {1,1, 1}. It should be 
noted, that the results obtained for this case are qualitatively 
identical to the more general case t = {£i,ti, t 2 } and v = 
{v 1 ,v 1 ,v 2 }. 

HI. SUPERLATTICE WITHOUT DISORDER 
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FIG. 1 : (Color online) Incompressible phases for an I = 3 superlat- 
tice with v = {y , ^, 0} and t — {1,1,1} without disorder. Solid 
line: iTEBD; crosses: DMRG; dashed line: CSCPE from @]. Simu- 
lation parameter for the iTEBD are x = 5, D - 1 = 3, fi/U = 1000 
(see appendix for definitions). DMRG results are obtained from a 
infinite size extrapolation. 



Let us consider first a superlattice BHM without disorder in 
the two cases v = {f,f,0}, t = {1,1,1}, and v = {0,0}, 
t = {1, 0.2}. As shown in Ref. H within a supercell mean- 
field approach such superlattices lead to loophole shaped in- 
sulator phases at fractional bosonic filling. In the following 
we will determine the boundaries of these loophole phases 
numerically and compare them to the mean-field predictions. 
Furthermore we will present a rather simple model which pro- 
vides an intuitive explanation. 

A. Numerical results 

In order to calculate the boundaries of the Mott phases for 
the superlattice Bose Hubbard model, we apply the iTEBD al- 
gorithm as described in the appendix to calculate the ground 
state of Hamiltonian ((TJ. We are able to calculate properties 
such as the local density g = (n) as a function of J and fi. 
Using this method to map out the shape of the insulating re- 
gions we make use of the fact that for any Mott phase the 
average local density is an exact multiple of l/l, I being the 
period of the superlattice. Thus the phase boundaries can be 
well approximated by the line where 

jttcr = where (n) = m/l ± e (3) 

for some m £ N (indicating the order of the lobe) and some 
l/l e > 1I22I1 . This line can be calculated by finding the 
value of fj, using a bisection method for a set of given m and 
J. Figures Q] and [2] show the results of this approach for the 
two different superlattice potentials specified above. 

The phase diagrams consist of a number of incompressible 
Mott phases, separated by a superfluid region. In contrast to 
the BHM for a simple lattice, there are however two types of 
insulating phases: The lobe-shaped ones, well known from the 
simple-lattice Bose-Hubbard model which have a finite extent 
at J = and the loophole-shaped ones which vanish at J = 0. 
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FIG. 2: (Color online) Incompressible phases for an I = 2 superlat- 
tice with v — {0, 0} and t — {1, 0.2} without disorder. Solid line: 
iTEBD; crosses with error bars: QMC (only for g = i) and dashed 
line: CSCPE (only for g = ±) both from GJ. Simulation parame- 
ter for the iTEBD are x ~ 7, D — 1 = 4 for g < 2 respectively 
D — 1 = 5 for g > 2, /3/U = 1000 (see appendix for definitions). 
For some phase boundaries e was set to 0.02 instead of 0.005 in order 
to avoid numerical artifacts. 

In general there are I distinct insulating regions for a superlat- 
tice of period A loophole is present, whenever the local 
potential vj is the same for two sites in the same unit cell. 
The following section will give a qualitative understanding of 
these loophole Mott regions. 

Figure Q] shows our results in the case v = {-j,-j,0}, 
t = {1,1,1}, together with the cell strong coupling per- 
tubative expansion (CSCPE) results from Q2D - As a check 
of our numerics we added numerical results from a density 
matrix renormalisation group (DMRG) calculation lfl3ll . 
The existence of non integer insulating phases at J = 
in this special superlattice has a direct connection to the 
case of a binary disorder-BHM 1 1 611 . which arises for exam- 
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pie in the presence of a second, immobile particle species 
(here of filling 1 with an inter species interaction of v = — ^). 

Figure [2] shows the numerical results together with the 
quantum Monte Carlo (QMC) results and the CSCPE data 
from U. The agreement between our numerics and the 
CSCPE is naturally good for small J but deteriorates for larger 
J. It is also apparent that while the insulator lobes are rather 
well described by the CSCPE approach, it is much less ac- 
curate for the loophole insulator regions, in particular for the 
case of varying potential depth (see figure[T|i. 



B. Two-site model 

We will argue in the following that the loophole insulator 
phases can entirely be understood from the effective dynamics 
within a unit cell of the superlattice. To this end let us discuss 
the above situation, where v = {i>i, vi, v 2 }, with v 2 < i>i. 

The presence of Mott lobes at fractional filling with a finite 
extend at J = can easily be understood along the lines of 
the simple-lattice BHM. As long as the filling is less than | 
the particles will occupy sites with local potential v 2 . Thus 
the chemical potential reads 



Mi = v 2 . 



(4) 



When the filling reaches the value i additional particles will 
start to occupy sites with local potential v\, giving rise to a 
particle-hole gap 



A/i 1 = fit - fii = vi 



V 2 . 



(5) 



To explain the existence of the loop-hole insulators, one has 
to take into account a finite hopping J. For J = any parti- 
cle added to the system between p = | and p = 1 increases 
the total energy by the same amount v-y . Thus the chemical 
potential stays the same. This picture changes however when 
a small but finite tunneling is included. If the filling exceeds 
the value | additional particles experience an effective super- 
lattice potential v = {vi,vi,U + v 2 } where the last term 
results from the interaction with particles already occupying 
sites with v 2 . If U > v\ — v 2 the superlattice effectively sepa- 
rates into degenerate double-well problems each correspond- 
ing to a unit cell. Due to the degeneracy of the double well any 
small tunneling J within the unit cell of the lattice needs to be 
taken into account while intra-cell tunneling can be ignored. 
A finite tunneling lifts the degeneracy of the single-particle 
states within the unit cell and leads to a splitting between 
symmetric and antisymmetric superpositions proportional to 
J. As long as the filling is less than p = | the particles 
occupy all sites with the smallest local potential v 2 and the 
symmetric superposition of the double-well {vi,vi}. After 
that additional particles have to go either to an already occu- 
pied side with potential v 2 , which is however suppressed by 
the large repulsive particle-particle interaction, or to the anti- 
symmetric superposition state. The latter requires an energy 
on the order of v\ + J, thus leading to another particle-hole 



gap on the order of J induced by intra-cell tunneling. More 
quantitatively the gap can be calculated by diagonalizing the 
two-site Hamiltonians H(Nb) for zero, one or two particles, 
(i.e. N B = 0, 1, 2), which read 



W(0) = 



in the basis {|00)}, 



H(l) 



V\ —J 
—J Vl 



(6) 



(7) 



in the basis {|10), |01)}, and 



H{2) 



U + 2V! -JV2 
-JV2 2«i -JV2 
-JV2 U + 2V! 



(8) 



in the basis {|20), |11) |02)}. The resulting ground state en- 
ergies are given by 



E(0) = 0, 
E(l) =v x -J, 



(9) 
(10) 



E{2) = -(u + Av 1 - Vl6J 2 + C/ 2 ) • (11) 

Calculating the chemical potentials fit = E(2) — £7(1) and 
fil = £7(1) - £7(0) yields: 



A*2 

3 

fit 



Vl — J 



= J - 

= Vl 



Vl ■ 

-J 



U 

~2 2 

0(J 2 ), 



u 2 



giving rise to a particle-hole gap 

Ap2 = J + 0(J 2 ). 



(12) 

(13) 
(14) 

(15) 



A generalisation of this discussion to the case of higher or- 
der loophole insulators or larger supercells is straight forward. 
For higher order loopholes, the accuracy becomes better, since 
the difference in the chemical potential between the rightmost 
and the other sites scales as U(n — 1), where n = [-yj + 1 is 
the number of particles of the corresponding Mott-insulating 
lobe. This means, that the effective two-site model gets even 
better for higher fillings. Qualitatively, one can understand the 
change in the shape of the loopholes for higher order (also see 
figure (3j just by considering the replacement J — > (n + 1) J 
in (O, because the single particle matrix H(l) is the only one 
relevant for the linear part of dT3b and because this replace- 
ment is the only influence of the other bosons already filling 
the lattice on the hopping in H(l). 



IV. SUPERLATTICE WITH DISORDER 

We now include a small disorder to the superlattice Bose 
Hubbard model. Of particular interest is the effect of the dis- 
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order to the loophole insulator phases. Disorder can be incor- 
porated into the model by replacing the last part of ([D accord- 
ing to: 



(16) 



with Aj being independent random numbers with continuous 
and bounded distribution, Aj £ [—A, A]. In the following 
we will restrict our analysis to the case of the superlattice 
potential v = {^,^'0}' * = { 1, 1, 1} and consider a 
canonical ensemble. 



A. Two-site model with disorder 

If the disorder is small, i.e. if 3A < U + 2u 2 — v\, the 
properties of the system in the vicinity of the A = loophole 
insulators can again be understood by considering the unit cell 
only, i.e. within an effective two-site model. The A = loop- 
hole insulator can be characterized by the number of particles 
per unit cell n and the disorder-modified chemical potential 
M = h'i + Ai,vi + A 2 }. 

Defining the total local energy in one unit cell as 



Tn u n 2 ■= ni(ni - 1) + ^n 2 {n 2 - 1) 



Ami 



the two-site Hamiltonians can be written as 

W n (2n-2) = T B _i, n _i J 
in the basis { \n — 1, n — 1 ) } for one particle less, 



n n (2n- 1) 



^n— l,n TlJ 
71 J ± n 71—1 



A 2 n 2 , 
(17) 

(18) 
(19) 



in the basis {\n — 1, n), |n, n — 1)} for zero extra particles , 
and 



n n (2n) 



(20) 







-Jy/n(n + l) 




T„. n -Jy/n(n + 1) 

-J^/n{n + 1) T„_i >n+ i 



in the basis {\n + l,n — 1), \n,n), \n — 1, n + 1}} for one 
additional particle. 

Calculating the chemical potentials one has to keep in mind, 
that the breaking up of the Mott-insulator is determined by the 
smallest particle-hole excitation throughout the whole system. 
Since all unit cells are decoupled, one has to find the disorder 
configuration which minimizes the energy gap. Therefore one 
has to calculate 



mm [^ 2 „(A 1 ,A 2 )-i; 2n _ 1 (A 1 ,A 2 ) 



(21) 



/il„_! = max [£ 2 „_i(Ai,A 2 )-£: 2 „_ 2 (Ai,A 2 )]. (22) 

3 Ai,A 2 



Since the disorder has only an effect on the local energy, we 
can set Ai = A 2 = A. The corresponding energies are given 
by 

£ 2 „_ 2 (A,A) = (n-l)(2A + C/-(n-2)), (23) 
S 2 „_i(A, A) = U(n - 1)2 + A(2n - 1) - nJ, (24) 
1 



£ 2n (A,A) = - [U {n - 1)2 + n{U + A A) 



-^&J 2 (n+l)n + U 2 \ 



(25) 



Minimization (maximization) of expressions (Bil l and d22b 
yields 



M 3n-1 
3 

Mt- 



-Jn + (n — 1)U + A, 



Jn 
1 



U 



Un 



y/8J 2 (n+l)n + U 2 . 



(26) 



(27) 



In order to calculate the critical tunneling rate at which the 
loophole insulator emerges, one needs to solve the equation 



/';' 



The solution can easily be found and reads 

U -2A 



Ji = A 



E/-4A 



(28) 



(29) 



for n = 1, and 



Jn = - 



1 



2n(n - 1) 



(U -AA)n 



-^n 2 (U 2 - 4UA + 8A 2 ) - 4A(U - 2A)r 



for n > 1. In both cases, the leading terms are given by 



Jn(A) 



iA + ^±±A 2 + 0(A 3 ), 



(30) 



(31) 



showing that the loophole decouples from the J = 0-axis in 
the presence of disorder, resulting in an insulating island. It 
should be noted that the two-site model cannot be used to cal- 
culate the maximum value of J for which the loophole insu- 
lator exists since for the vanishing of the gap at large J values 
also inter-cell tunneling processes need to be taken into ac- 
count. 



B. Numerical results 

As seen above from the effective two-site model, the loop- 
hole insulator regions are decoupled from the J = axis 
giving incompressible islands. Although the system is for 
any given disorder realization not translational invariant, the 
iTEBD method can be used also in this case. To this end we 
define supercells each of which with the same disorder. The 
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FIG. 3: (Color online) iTEBD results for boundaries of incom- 
pressible phases for I = 3 superlattice with v = {tt, y,0} and 
t = {1,1,1} and a disorder amplitude A/ (7 = 0.04. Dashed line: 
pure case, solid line: disordered case. For simulation parameters see 
figure [TJ Horizontal lines indicate the positions of density cuts of 
figure[4] 



supercells have to be large enough such that effects from spa- 
tial correlations and finite size can be ignored. We have cho- 
sen a supercell length of 96. Increasing this length did not 
show any noticeable changes. To calculate the physical quan- 
tities, the numerical results have to be averaged over a number 
of different disorder realisations, namely over different sets of 
disorder A = {Ai, A 2 , . . . , A/} with Aj £ [-A, A], where 
a boxed disorder distribution is assumed. The length of the 
vector A is the same as the size of the system simulated, see 
appendix. It turns out that 20 realizations provide sufficient 
convergence for the purpose of this paper. 

Figure [3] shows the results of the iTEBD calculations both 
in the pure (A = 0) and the disordered (A = 0.04C/) case. 
The first thing to notice is the shrinking of the Mott-insulating 
lobes for J = due to the disorder. As known from the BHM 
ifll [l5tl the Mott-lobes shrink by an amount of 2A at the 
J = axis. The second and more important thing to notice is 
the decoupling of the loophole insulator from the J = axis, 
meaning that there is no insulating phase for the respective 
filling for J < J cr it, in full agreement with equations (|29l[30l l. 

The decoupling of the loophole insulator from the J = 
axis can most easily be seen in a cut parallel to the /i-axis for 
fixed J, showing the average local density as a function of 
the chemical potential. In the case of small hopping without 
disorder, this cut shows, beside the expected Mott-lobes 
at g = | and g = 1, an intermediate plateau at filling |, 
corresponding to the loophole phase (Figure [4] lower plot, 
solid line). For larger hopping (upper plot, solid line), the 
width of the plateau is slightly increased according to the 
shape of the loophole in figure Q] In the case of disorder, the 
plateau for g = | vanishes for small hopping (figure|4] lower 
plot, dashed line). For large hopping (upper plot, dashed line), 
the incompressible phase survives, however with a largely 
reduced width compared to the pure case, which shows the 
decoupling of the loophole from the J = 0-axis as predicted. 
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FIG. 4: (Color online) Density cut along the horizontal lines in figure 
[3] for the pure (solid line) and the disordered (A/U = 0.04, dashed 
line) case. Upper plot: cut along the upper line in figure[3]at J/U = 
0.06, lower plot: cut along the lower line in figure|3]at J/U = 0.02. 
Simulation parameters see figureQ] 



In figure [5] we show the numerical results for the first, sec- 
ond and third loophole for increasing values of the normalised 
disorder amplitude. It should be noted that due to the finite 
number of disorder realisations it is difficult to accurately de- 
termine the lower tip of the insulating island in the numerics. 
In figure [6] we compare the onset of the insulating loopholes 
obtained from the analytic two-site model with numerical re- 
sults. One immediately recognizes two things: First, the onset 
of the loophole, J cr it, is a monotonous function of A and sec- 
ond, the higher the filling of the lobe, the earlier the insulating 
region arises. The numerical value of J cr i t was obtained by 
reading off the values from the numerically determined phase 
diagram assuming generous error margins l23ll . Taking into 
account these errors, figure [6] shows a rather good agreement 
of the two-site model to the numerics. However, the analytic 
prediction tend to be slightly to large compared to the numer- 
ics, nevertheless giving the right leading order for small A. 
This is because for larger A the critical hopping gets larger 
than allowed by the assumption of a decoupled two-site prob- 
lem. By diagonalising the complete three-site unit cell with 
periodic boundary-conditions, which gives another but less in- 
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FIG. 5: (Color online) Detailed analysis of the first three loophole 
insulators (from left to right: g — 2/3, q = 5/3, q = 7/3) with 
varying disorder amplitude, increasing from the outer to the inner 
lines (black: A/U = 0.00, magenta: A/U = 0.02, red: A/U = 
0.04, orange: A/U = 0.06). For the simulation parameters see 
figure[T]except for D — 1 = 4 in the rightmost plot (q = 8/3). 
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FIG. 6: (Color online) Critical point J n of the onset of the loophole 
insulating island as a function of the disorder. Solid line: analytic 
prediction from equations i29\ and i BOI l. crosses: data with error bars 
as read from figure[5](not all used data is shown in figure|5](, dashed 
lines: leading order in eq. d31b . From top to bottom: n = 1, n = 2, 
n — 3. 



ular the case of degenerate potential energies and/or degen- 
erate tunneling rates within the unit cell of the superlattice 
have been discussed. Using both, exact numerical meth- 
ods such as the infinite-size time evolving block decima- 
tion (iTEBD) algorithm, and the density matrix renormaliza- 
tion group (DMRG) we calculated the boundaries of incom- 
pressible Mott-insulating phases. The existence of additional 
loophole-shaped Mott domains, predicted before, was veri- 
fied and their numerically determined phase boundaries com- 
pared to other approaches such as the cell strong coupling ex- 
pansion. A simple effective model was presented that takes 
the full dynamics within a unit cell into account. The model 
provides a rather straight forward explanation for the emer- 
gence of loophole Mott domains in the case without disorder. 
Adding a small amount of disorder with continuous, bounded 
distribution lead to a shrinking of the loopholes to Mott insu- 
lating islands with the remarkable feature of a compressible 
to insulating transition with increasing bosonic hopping. The 
analytic predictions for the critical hopping for this transition 
from the effective model were compared to numerical simula- 
tions and found in very good agreement. 
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APPENDIX: THE TEBD ALGORITHM AND THE ITEBD 
IDEA 



tuitive approximation, we get a curve for J cr it that is below the 
two-site prediction, but is the same in first order and in better 
agreement with our numerics. 

V. SUMMARY 

In the present paper we have analyzed the superlattice 
Bose-Hubbard model with and without disorder. In partic- 



In the following we give a short summary of the numeri- 
cal algorithm used in sections fill Al andl lV Bl The basic idea 
of the TEBD algorithm emerged from quantum information 
theory fill] and it can be used to simulate one-dimensional 
quantum computations that involve only a limited amount of 
entanglement. 

Here we want to use it for an imaginary time evolution of 
the one dimensional Bose Hubbard model. The state of the 
system can be represented as a matrix product state: 



X D 

m= V V r [1]ii A [11 r [21i2 x lL - sl T [L - 1]iL - 1 \ [L - 1] rW iL h, u) (32) 

l w / 2-~i L^t "i Q i a i Q 2 ■ ■ ■ A a L - 2 1 at L -20tL-i V^i 1 ctL-i I* 1 ■ ■ ■ lL >- yJZ -> 

a 1 ,a 2 ,...a i =l i 1 ,i2,...i L =0 
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Here D is the dimension of the local Hilbert space on a sin- 
gle site and \ is the number of basis states in the Schmidt 
decompositions (see below) to be taken into account. \ is a 
measure for the maximum entanglement in the system and is 
assumed not to increase with system size or to increase only 
very slowly. In our model the number of particles per site is in 
principle not bounded. But to reduce the numerical effort one 
can safely set a maximum number of particles D — 1 allowed 
per site, since higher occupancies are strongly suppressed due 
to the on-site interaction. \ii . . . %l) is the state from the Fock 
basis, where there are ik bosons on site k. 

Furthermore we require our matrix product representation 
to be be in the canonical form, i.e. equation (f32t represents 
the Schmidt decomposition for any bipartite splitting of the 
system at the same time. This means for any given k, the 
Schmidt decomposition between sites k and k + 1 is given by 



I*) = 



E A L fe W 

a=l 



(33) 



where the Schmidt coefficients Aq ' are normalised as 



(34) 



andthe {|*£ " fe ~ 1] }} ({ |*L fe "' L] )j ) form an orthonor- 

mal set of states in the subspace of the first k (last L — k) 
sites. By sorting the Schmidt coefficients in an non ascending 



order for every bond, this makes the representation de facto 
unique. Explicitly, 



E 



r[i]«i 



\[k-l] T [k] lk I- ■ \ 



(35) 

where the A's account for all the Schmidt coefficient and the 
A's care for the transformation into Fock space at every sin- 
gle site. Describing an arbitrary state in general requires that 
the Schmidt number \ is of the order D L . We will use how- 
ever a relatively small, constant \ to avoid exponentially in- 
creasing c omp lexity of the numerical problem. It has been 
shown in 11711 . that this seemingly strong assumption is jus- 
tified and gives a good approximation for the ground state. 
The latter is related to the fact that the ground state of one- 
dimensional systems with finite-range interactions has either 
a constant entanglement (for noncritical systems) or the en- 
tanglement increases only logarithmically with the size (for 
critical systems). Small values of \ give usually very good re- 
sults for local observables, while correlations are only poorly 
approximated over very large distances. For the latter the ap- 
proximation can be improved by choosing a larger \ propor- 
tional to the distance |12|1 . The amount of coefficients needed 
to specify the matrix product state with given fixed \ is of the 
order L ■ D ■ \ 2 an d can be handled numerically in contrast 
to the D L coefficients required for representation in the full 
Fock space. 

Expressing the state in a local basis for sites k and k + 1 
only, 



i*>= E E^ llr W lr ^ lw ^ 1] i*^ ,A " ll >wiJ>i*^-" LI ) 

ck,/3,7=1 i,j—0 



(36) 



we see that applying an operator that involves sites k and k + 1 
only is equivalent to manipulating the matrices r[ fc ] and r[ fe+1 l 
and the vector only which is implemented as follows. 
For reasons of stability we use A^g := r^l 1 Ao* through- 
out the algorithm ll24ll . To shorten the notation we rename 



aL*- 11 |*£- 

giving 



.fe-H 



and A 



[fe+l]|^[fc+2...Z,h 



I7) in (|36]l 



|*> = EE^ 



[k]i .[k+m 

a/3 



a 7 
ij 



[k+1 



■\aijj). 



(37) 



The objective is now to decompose T into a product of 
matrices A^Ap^ 1 ^ and to keep the canonical form. The 

are the eigenvectors of the reduced density ma- 







trix 



j[k+i-L] _ Xr [1 '" fc ]|^)(*| 



EE^ 11 



071 



3132 «* 
7172 v 



M J132 
m 7172 



rplj 2 ] 
a~f 2 ) 



(39) 



b'i7i) O272 1 
,[fe+l] x [fc+l] 

A 7l A ~f2 



Applying a two site operator V given by the matrix V^ n then 
results in 



«7 Imp 



(38) 



Diagonalising M gives the new A 



the new ( Ag ' ) as eigenvalues. (Using the T matrices in- 



[k+l]j 



as eigenvectors and 



stead of the A matrices would require a division by A 



[k+i] 



But A 



[k+i] 



can be zero if the Schmidt number for this bond 



is smaller than In general there are D ■ \ nonzero eigen- 
values. (This is due to the possible creation of entanglement 



8 



by V.) But we can only keep the x biggest of them. There- 
for we have to renormalise the new A^ according to fl34l l. 
This is necessary anyway if we have a non-unitary V as in the 
case of an imaginary time evolution. The are given by 

In order to calculate the ground state of our system (see 
lfl9l0 . we divide the Hamiltonian (fl3 into two parts H CV cn and 
H dd, where H cvcn (H dd) couples sites j and j + 1 for even 
(odd) j only. The local parts of H can be distributed between 
H cvcn and H D dd arbitrarily. The ground state is then given by 
an imaginary time evolution 



I Aground) = lim ■ 
p— *oo 



-^1*0) 



(40) 



Here any initial state |^o) is sufficient, as long as it has a finite 
overlap with the (yet unknown) ground state. The evolution is 
implemented by repeatedly applying small time steps er He ', 
so called Trotter steps. The norm is conserved in this pro- 
cedure (see above). So after T steps only the ground state 
has a reasonable contribution to our state \f (3 = T ■ e is 
much bigger than the inverse of the energy of the first ex- 
cited state (relative to the ground state energy). In order to 
write e~ He as a product of two site operators we use the 



Suzuki-Trotter decomposition 11811 . In first order one can get 

e -He _ e -H omn e e -H odd e _|_ 0(e 2 ), in second order e~ He = 

e s f EiLe e -i ^ odde e W 2 - 6 + 0(e 3 ). For higher orders see 
lfl8ll . Thus we can calculate the ground state by repeatedly 
applying two-site operators. 

To calculate expectation values of observables we again 
take a look at d36b . The expectation value of a nearest neigh- 
bour observable, say (^|ata.fc+i|\I r ) can be directly calculated 
because all occurring states are mutually orthogonal and nor- 
malised. An nth site nearest neighbour observable can be cal- 
culated by expressing the state in the local basis for site k to 
k + n analogous to (l36l l. For non nearest neighbour observ- 
ables we can use the swap gate to bring the sites of interest 
together ifTHl. 

A powerful feature of the algorithm is its application to infi- 
nite, translationally invariant systems. Suppose a Hamiltonian 
that has a periodicity of c sites (as (Q~|i has for a superlattice), 
restricting to c = 2 for clarity. The state of an infinite system 
is a slight modification of ( 1321 . 



I*) 



\" \[k-l] r [k]i k A [fc] r [fc+l]i fc+ i A [fc+l] i- 
"fc-i &k-l<*k <*k a k "afc+i ■■•1"- 



(41) 

The imaginary time evolution is started with a translationally 
invariant state, so all T's and A's are the same in the begin- 
ning. The scheme in figure [7] shows, that the c-periodicity of 
the representation is preserved during real or imaginary time 
evolution. This is because all two-site operations H dd and 
Hcvcn are the same respectively and are all applied to every 
other pair of matrices. 
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FIG. 7: (Color online) Symbolic representation of the effect of the 
TEBD algorithm to a translationally invariant state. The uppermost 
(initial) state is first changed by application of H a dd, giving new T's 
and A's. The second step with H even then produces a list of alternat- 
ing T's and A's such that only two of them need to be kept in memory. 
Every further Trotter step preserves this symmetry. 



So we only have to store two T matrices and two A vectors. 
It is even more important that we only have to apply two two- 
site operators per Trotter step. 

After imaginary time evolution we end up with an c- 
periodic ground state. (That means that expectation values 
have a periodicity of c sites. Although there can be contri- 
butions in fiTT i from states which have a nonperiodic Fock 
representation. This is a clear distinction from the case of 
periodic boundary conditions, where not only all expectation 
values, but also the wave function must be periodic.) This 
is called the iTEBD algorithm Il2ll . This means that we can 
efficiently calculate observables in the thermodynamic limit. 
If we were using DMRG or normal TEBD we would have to 
simulate large finite systems, which is time consuming, and 
then extrapolate to L = oo to get rid of finite size effects but 
introducing additional error. 

The idea works as well for c > 2. If c is odd, we have to 
choose 2c as period, since we need a clear distinction between 
Hcvcn and H Q dd- In fact we used it in this work for the non pe- 
riodic Hamiltonian of the disordered superlattice model, thus 
not saving calculation time (a large value has to be used for c 
in order to have a sufficiently random disorder) but getting rid 
of boundary effects. 

^ Finally we note that the TEBD algorithm itself is in prin- 
ciple only correct for unitary operations. Non unitary oper- 
ations were found to destroy the representation in the sense 
that the Schmidt vectors in d3~3l are no longer exactly orthog- 
onal, i.e. the representation is no longer canonical 12011 . Ad- 
ditional steps to conserve orthogonality in the algorithm were 
proposed in l2lll . These were not incorporated here, since for 
small e the Trotter steps are quasi orthogonal. Numerical anal- 
ysis shows, that the scalar products of the normalised Schmidt 
vectors in the resulting ground state are of the order 10 -3 . 



9 



[1] D. Jaksch, C. Bruder, J. I. Cirac, C. W. Gardiner, and P. Zoller, 

Phys. Rev. Lett., 81, 3108 (1998). 
[2] M. Greiner, O. Mandel, T. Esslinger, T.W. Hansch, and I. Bloch, 

Nature 415, 39 (2002). 
[3] L. Santos, M. A. Baranov, J. I. Cirac, H.-U. Everts, H. 

Fehrmann and M. Lewenstein, Phys. Rev. Lett., 93, 030601 

(2004) 

[4] S. Peil, J. V. Porto, B. Laburthe Tolra, J. M. Obrecht, B. E. 

King, M. Subbotin, S. L. Rolston, and W. D. Phillips, Phys. 

Rev. A 67, 051603(R) (2003). 
[5] V. G. Rousseau, D. P. Arovas, M. Rigol, F. Hebert, G. G. 

Batrouni, and R. T. Scalettar, Phys. Rev. B 73, 174516 (2006) 
[6] R. Roth, and K. Burnett, Phys. Rev. A, 68, 023604 (2003). 
[7] P. Buonsante and A. Vezzani, Phys. Rev. A 70, 033608(R) 

(2004). 

[8] P. Buonsante, V. Penna and A. Vezzani, Phys. Rev. A 70, 

061603(R) (2004). 
[9] P. Buonsante and A. Vezzani, Phys. Rev. A 72, 013614 (2005). 
[10] L. Fallani, J. E. Lye, V. Guarrera, C. Fort, and M. Inguscio, 

Phys. Rev. Lett. 98, 130404 (2007). 



[11] G. Vidal, Phys. Rev. Lett. 91, 147902 (2003). 

[12] G. Vidal, Phys. Rev. Lett. 98, 070201 (2007). 

[13] U. Schollwock, Rev. Mod. Phys. 77, 000259 (2005). 

[14] M. P. A. Fisher, P. B. Weichman, G. Grinstein, and D. S. Fisher, 

Phys. Rev. B 40, 546 (1989). 
[15] J. K. Freericks and H. Monien, Phys. Rev. B 53, 2691 (1996). 
[16] A. Mering and M. Fleischhauer, Phys. Rev. A, 77, 023601 

(2008). 

[17] F. Verstraete and J. I. Cirac, Phys. Rev. B 73, 094423 (2006). 
[18] M. Suzuki, Phys. Lett. A, 146, 319 (1990). 
[19] G. Vidal, Phys. Re v. Lett. 93, 040502 (2004). 
[20] R. Orus, G. Vidal, |arXiv:0711.3960l (2007). 
[21] Y.-Y. Shi, L.-M. Duan, G. Vidal, Phys. Rev. A 74, 022320 
(2006). 

[22] Throughout the paper we set e = 0.005 

[23] Calculating the tip of a lobe is in general numerically expensive. 
[24] The vectors A' fe ' have all to be kept separately in order to not 
loose information about the canonical form of the state. 



